#ifndef _BEM_UTIL_H_
#define _BEM_UTIL_H_

#include <mat.h>
#include <vec.h>

Mat22d evaluateGradientQuadraturePoint(Vec2d p, Vec2d q, Vec2d qn);
Mat22d evaluateQuadraturePoint(Vec2d p, Vec2d q);

double tensorTensorContraction(Mat22d A, Mat22d B);
Mat22d tensorVectorContraction(double A[2][2][2], Vec2d v);
Mat22d tensorVectorDot(double A[2][2][2], Vec2d v);
double* tensorProductGradient(double s, Vec2d ds, Mat22d A, double dA[2][2][2]);
double* getDrI(Vec2d r, double rn, Mat22d I);
double* getDrM(Vec2d r, double rn, Mat22d M);

double* getKernelGradient(Vec2d p, Vec2d q);

int getIndex(int i, int j, int k);
int getTensorArrayIndex(int i, int j, int k);

#endif